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Abstract 

We present a study of the effect of artificial 
dissipation models on nonlinear wave computations 
using a few high order schemes. Our motivation 
is to assess the effectiveness of artificial dissipation 
models for their suitability for aeroacoustic compu- 
tations. We solve three model problems in one di- 
mension using the Euler equations. Initial condi- 
tions are chosen to generate nonlinear waves in the 
computational domain. We examine various dissi- 
pation models in central difference schemes such as 
the Dispersion Relation Preserving (DRP) scheme 
and the standard fourth and sixth order schemes. 
We also make a similar study with the fourth order 
MacCormack scheme due to Gottieb and Turkel. 

KEY WORDS: Artificial dissipation, High or- 
der schemes, Computational aeroacoustics 

1. Introduction 

Very high numerical accuracy is essential in 
aeroacoustics computations. To get such accuracy, 
one should use high order schemes with very lit- 
tle numerical dispersion and dissipation. Numeri- 
cal solutions with high order schemes are usually 
very good in flows without shocks. However, be- 
cause of inadequate numerical dissipation a high 
order scheme may generate spurious oscillations in 
presence of shocks or steep gradients. In this study 
we consider computations of flows with very sharp 
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gradients. Such flow conditions may appear in su- 
personic jet noise computations. All basic schemes 
used in this study are formally at least fourth order 
accurate in space. We augment these basic schemes 
with various artificial dissipation models and eval- 
uate these artificial dissipation models. Tam and 
his colleagues specifically designed basic scheme for 
aeroacoustic calculations and it is known as the Dis- 
persion Relation Preserving (DRP) scheme*’ 3 
The DRP is a central finite difference scheme. We 
also consider the standard fourth and the sixth order 
central difference schemes and the fourth order Mac- 
Cormack scheme due to Gottlieb and Turkel"*. The 
fourth order MacCormack scheme has been tested 

r 

for accuracy in aeroacoustic applications 0 ’ and is 
extensively used, e.g.^’ *^ . 

In their study of the DRP schemes for non- 
linear wave computations, Tam and his colleagues 
discovered that the wavelengths of spurious oscil- 
lations are concentrated in a narrow band in short 
wave range. They formulated a selective artificial 
damping model to suppress such oscillations and 
used it for computation of nonlinear pulses^ . There 
are many other dissipation models which are also 
shown to be effective in flow computations. One 
such dissipation model was proposed by J ameson et 
al. **. This model is widely used in aerodynamic 
calculations. Turkel* ^ extended this model to a set 
of matrix- valued coefficients to give appropriate vis- 
cosity for each wave component. In this study, we 
will consider these dissipation models and also the 
model proposed by MacCormack and Baldwin 1 . 

In the next section we give the governing equa- 
tions. Descriptions of the basic schemes and arti- 
ficial dissipation models are presented in sections 3 
and 4 respectively. We give our test problems in 
section 5 and discuss results in section 6. 
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2. Governing Equations 

We solve the one dimensional Euler equations 
written in the following form 

Qt + fx — 0 

where 



the Laplace transform of the finite difference scheme 
is a good approximation of the partial derivative. 
Their scheme can be written as 

= --k £ * (» i«) 

>==-3 

3 

Q? +1 = Q: + At bjK*-’ (3.16) 

i=o 


/-LT + , 

\ puH ) 

P = (V~ 1 )(E ~ \( m2 ) 

P 

where p, u, p, 7, E and H are the density, velocity, 
pressure, ratio of specific heats, total energy and 
enthalpy respectively. 

3. Basic Schemes 


oq = 0; a! = — a_i = 0.79926643 

a 2 = — a _2 = -0.18941314 

03 — — a _3 — 0.02651995 

6 0 = 2.30255809; 61 = -2.49100760 

&2 = 1.57434093; 63 = -0.38589142 

Above DRP scheme is formally fourth order accu- 
rate in space and third order accurate in time. In 
a later study, Tam and Shen^ revised the values of 
the coefficients aj to obtain a better overall numer- 
ical wave characteristics for a seven point stencil. 
Revised aj coefficients are 


We consider a high order MacCormack scheme 
along with a few central difference schemes. The 
MacCormack scheme uses forward and backward 
differences in predictor and corrector sweeps. All 
basic schemes presented in this study can be writ- 
ten in conservative discrete form. Such conserva- 
tion is very important for flows with shocks and is 
illustrated by the fundamental theorem of Lax and 
Wendroff^ * 


3.1 The DRP and Central Schemes 

Tam and Webb* developed a class of high or- 
der scheme known as the Dispersion Relation Pre- 
serving (DRP) scheme for computational acoustics. 
They used centred differencing for spatial discretiza- 
tion. The coefficients of their spatial discretization 
scheme were chosen by requiring the Fourier trans- 
form of the finite difference scheme be a close ap- 
proximation to that of the partial derivatives. Tam 
and Webb used a seven point stencil for fourth order 
spatial accuracy. They proposed a multilevel finite 
difference time integration scheme. Coefficients of 
this time integration algorithm were chosen so that 


ai = — a_i = 0.770882380518 
a 2 = — a _2 = -0.166705904415 
a s = -a _ 8 = 0.020843142770 

For all computations with the DRP scheme in 
this paper, we used the above set of values for aj. 
This DRP scheme has a seven point stencil and for- 
mal fourth order spatial accuracy. The coefficients 
aj in the standard fourth order central difference 
scheme are 


*o = 0 ; 



— a_ 2 = — 


1 
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The stencil is five points wide, i.e., as — a_s= 0. 
Similarly the standard sixth order central difference 
scheme on a seven point stencil is 


ao = 0; 



*2 


— a_2 = — 


3 

20 


as = — a_ s = 


60 
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For time integration of the DRP, the standard 
fourth and the sixth order central difference schemes 
we use the algorithm given in eq. 3.1b. Thus the 
formal temporal accuracy for these calculations is 
third order. 

3.2 High order MacCormack Scheme 


predictor and corrector steps as in the case of Mac- 
Cormack Baldwin model (§4.2). 

4. Artificial Dissipation Models 
4.1 Tam’s Model 


Gottlieb and Turkel^ extended the standard 
MacCormack scheme (2nd order accurate in both 
space and time) to a spatially fourth order accurate 
scheme. In this study we will refer this scheme as 
the fourth order MacCormack scheme. This scheme 
has a predictor and a corrector stage and for one 
dimensional computations may be written as: 

the predictor step with forward differences is 


<?; = <?? + 6^{7(/T + 1 - f?) - (/? +2 - /T+i)} 

(3.2a); 

and the corrector step with backward differences is 

Q»” +1 = \[Qi + <?? 

+ g^{7 (fi - fi-i) - (fi-i - fi- 2 )}] (3-26) 

This scheme is second order accurate in time and 
becomes fourth order accurate in the spatial deriva- 
tives when alternated with symmetric variants 
We define Li as a one dimensional operator with 
a forward difference in the predictor and a back- 
ward difference in corrector. Its symmetric variant 
L 2 uses a backward difference in the predictor and 
a forward difference in the corrector. Therefore to 
ensure the fourth order spatial accuracy, the sweeps 
are arranged in our computations as 

Q n+1 = LiQ n 
= L 2 Q n+1 

In this study when we add any disspation model 
in this scheme (except the MacCormack- Baldwin 
model in §4.2) we modify only the corrector step 
as 



Q” +1 = \[Qi + Qi 

+ e ^ {7(/< ~ - 1*- 1 - *- 2 » ] + AtD ' 

where Di is the artificial dissipation at location i. 
Alternatively, one can also add dissipation in both 


The basic DRP scheme (§3.1) has very little 
numerical disspation. In order to compute nonlinear 
waves with this scheme, Tam and his colleagues^ > ^ 
formulated an artificial selective damping term to 
be added to the basic algorithm. Thus, the eq. 3.1a 
becomes 


*"=-xr £ <■>*> + 


Ax / o 

j = — 3 


D? = - 


AXR ‘ > = — 3 


E 


+jf 


(4.1) 


where u * = — ti™” | is the difference between 

the maximum and the minimum velocity in the com- 
putational stencil centered at location t. The coef- 
ficients Cj are determined by requiring the damping 
function be a close approximation to a Gaussian 
function. Tam and Shen^ recommended following 
numerical parameters after extensive numerical ex- 
periments: 


c 0 = 0.327698660846; c x = -0.235718815308 

c 2 = 0.086150669577; c 3 = -0.014281184692 

The damping function (D) is an even function, i.e., 
c_j = Cj. Their recommended value for the sten- 
cil Reynolds number (£,) is 0.1. Recently Tam^ 
used the stencil Reynolds number equal to 0.05 in 
his computation of the wave equation. We would 
like to point out that Tam and his colleagues for- 
mulated this dissipation model to specifically main- 
tain the correct propagation characteristics of waves 
and damp short waves generated by numerical so- 
lution, rather than to maintain high order of accu- 
racy in the sense of expanded Taylor series. Thus 
this model should yield good results related to wave 
propagation characteristics. 


4.2 MacCormack- Baldwin Model 

MacCormack and Baldwin^. ^ devised a dis- 
sipation model for the MacCormack scheme. Thus, 
the modified flux becomes 
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4.4 Turkel’s Matrix- Valued Model 


Si 


= /« + e(H + 


c) |Pi+i ~ 2 P» .+.Pl-i j d Q^ 
Pi+i + 2 &+ jh. i 


where dQi is a forward difference (i.e., Qi+ 1 - Qi) 
in the backward sweep and vice versa . The physical 
flux and the speed of sound are / and c respectively. 
The modified scheme with artificial viscosity is ob- 
tained by replacing / by / in the basic scheme (i.e. 
in eq. 3.2a and 3.2b for the fourth order MacCor- 
mack scheme). The parameter e should be between 
0 and | for stability. In our tests we use e = 


Turkd^ proposed a matrix scaling instead of 
scalar scaling in the JST dissipation model. This 
matrix dissipation model was tested by Jorgenson 
and Turkel^. One major objective of this model 
is to reduce the amount of artificial dissipation for 
slower waves. Thus the scaling parameter A in eq. 
4.3 becomes a matrix. Details of this matrix are 
given in the above cited references. We implemented 
this model as follows: 

Ax 


4.3 JST Model 

Jameson, Schmidt and Turkel* 1 proposed an 
artificial dissipation model where the convective flux 
is modified as 

/«+! = fi+\ ~ 

where 

d t+ > = \ + >[eQdW t - e W(dW t+1 - 2dW t +dW i _ l )] 

(4.3) 

dWi = W i+1 - Wi 

= * (2) max(i/ i+1 , i/<) 

= max(0j*< 4 > -£>,]) 

_ |Pi+i - 4-p»-i| 

J*+i + 2pi +pi_i 

w = Q + [0,0,p] T 

kS 2 ) = i, /c( 4 ) = and A is a scaling factoi. In 
this study we use A = |ti| + c. The second differ- 
ence term adds dissipation near the shock while the 
fourth difference term damps high-frequency modes 
and reduces to zero near the shock. We will refer 
this model as the JST dissipation model in this pa- 
per. The damping function (A) for this model is 
as follows: 


A = \ 3 I+{^±±1 - As )[ + G 2 ] 
l c 



Rx = (y,-«,l) 
— ( _ lj 0) 

Rs = (« 2 , l) 




|tt+ c\ 1 

A 2 

— 

|«-c| mai(n, n+i) 

k As , 


. M max(/tj,^ + i) , 


Pi+1 +2p< +Pj_i. 


Di = 



Ax 


|T <+1 -2T < + T < - 1 | 

T i+ 1 + 2 Ti + Ti_! 
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dQ, — <?i+i Qi 
e W(j) = K W\ j 

eW(j) = max[0, - e (3 >(j)] 

AQi(j) = «r ) (i)dQ i +el 4) (i)(dQ i+ i-2dQ i +dQ i _ 1 ) 

where T is temperature, parameters k W and 
are same as in the JST model. 

5, Test Problems 

We consider three test problems to evaluate 
various dissipation models in the basic schemes. 
The first two problems were previously used to eval- 
uate various algorithms in the ICASE/LaRC Work- 
shop on Benchmark Problems in Computational 
Aeroacoutic^ . In addition we consider Lax’s shock 
tube problem. 

5-1 Test Problem 1 

In this problem a initial pulse steepens with 
time. Initial conditions are given as: 

U = i e - (*) a ( ,n2 ) 

2 



P = ( 1 + 

The computational domain is -50 < x < 350. 


5-2 Test Problem 2 

This is similar to a standard shock tube prob- 
lem. There is a very sharp gradient instead of a 
discontinuity in the initial condition. Initial condi- 
tions are 


u 

= 0 


p 

= 4.4 

x < -2 

V = 

2.7 + 1.7co*[te±2Jl], 

-2 < x < 2 

P 

= 1 . 

x > 2 


p = (tp) ’ 

The computational domain is -100 < x < 100. 


5.3 Test Problem 3 (Lax’s problem) 

This is a Riemann (shock tube) problem used 
by Lax^. Initial conditions are 

u = .698, p = 3.528, p = .445 for 0 < x < .5 
u = 0, p = .571, p = .5 for .5 < * < 1 
The computational domain is 0< x < 1. 

6. Results 

One uses artificial dissipation to eliminate spu- 
rious numerical oscillations in computed flows with 
shocks or steep gradients. Too little dissipation may 
not be enough to remove spurious oscillations, while 
too much dissipation will smear of the solution pro- 
file and destroy the accuracy of the computed so- 
lution. Therefore one should try to add appropri- 
ate amount of dissipation to a basic scheme. Our 
primary goal in this study is to examine the effec- 
tiveness of a few artificial dissipation models. We 
tried to keep the model constants same as in cited 
references and did not attempt to optimize any con- 
stant. We tested these dissipation models in a few 
basic schemes. A comparison of two such schemes, 
namely the DRP and the fourth order MacCormack 
scheme (without dissipation) for the first two test 
problems was presented in an earlier study^. Here 
we also consider the standard fourth and sixth or- 
der central difference schemes. For all our com- 
putations reported in this paper we used 401 grid 
points and the characteristic boundary condition at 
the boundaries. In the first test problem, an ini- 
tial pulse steepens with time and nonlinear effects 
become very important. In Figure 1, we show the 
evolution of the density profile for this test prob- 
lem. These computations were made using the DRP 
scheme with the Tam dissipation model with the 
stencil Reynolds number (R,) equal to 0.05. We 
observe that the initial profile steepens quickly and 
at t = 100 the leading edge becomes very sharp. We 
chose this time level to evaluate the damping models 
and the basic schemes for computations of the test 
problem 1. All computations with the central differ- 
ence schemes (the DRP, standard 4th and 6th order 
central difference scheme) presented in this paper 
were done with CFL number equal to A_. In Fig- 
ure 2, we show density profiles computed using the 
DRP scheme at t=100. There is excessive amount 
of oscillations in the computed solution when no ar- 
tificial dissipation is used. Both the Tam dissipa- 
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tion model (i?, = 0.1) and Turkel’s matrix valued 
disspation model eliminated most of the oscillations. 
We used «( 2 ) = 4 and k W = ^ for computations 
with the matrix valued dissipation model. These 
values of #c are 16 times the recommended values 
in the JST model (§4.3). As indicated earlier, par- 
ticular values of the constants may be optimized 
for a given basic scheme. Here we chose the values 
of k only to ensure sufficient damping. Parametric 
studies of the constants may yield a better choice. 
In addition to different choices of one may also 
try to optimize this dissipation model with differ- 
ent choices of switching parameters e.g., Vi, etc. 
in §4.4. However, fine tuning of the parameters is 
expensive and may not be always feasible in practi- 
cal computations. The Tam dissipation model has 
a scaling parameter known as the stencil Reynolds 
number(i2«). Effect of this parameter and also a 
comparison of the standard 4th and the 6th order 
central difference schemes with the DRP scheme are 
shown in Figure 3. The DRP scheme which we con- 
sider in our tests is formally 4th order accurate in 
space and uses a seven point stencil. The standard 
4th order central difference scheme uses a five point 
stencil and the standard 6th order central difference 
scheme uses a seven point stencil. We observe that 
the DRP scheme performs better than the standard 
4th order scheme. However, the differences between 
the DRP scheme and the standard 6th order central 
difference scheme for our present computations ap- 
pear to be very small. These results show the same 
trend as in wave equation computations presented 
by Tam^. The artificial dissipation may have re- 
duced the differences among the schemes. We also 
observe that with the change of the stencil Reynolds 
number (R,), both the solution profile and the shock 
location change slightly. 

Next we evaluate the dissipation models for 
computation of the test problem 1 with the fourth 
order MacCormack scheme. We used the CFL num- 
ber equal to | for all computations in this paper 
using the fourth order MacCormack scheme. Den- 
sity profile computed using this scheme without any 
added artificial dissipation is shown in Figure 4. 
We also show the computed profile using the DRP 
scheme augmented with the Tam dissipation model 
( R t = 0.1). The fourth order MacCormack scheme 
with its builtin dissipation appears to generate mi- 
nor spurious oscillations. One may use artificial dis- 
sipation model to reduce these oscillations. Effects 
of various dissipation models on the computed den- 
sity profiles near the shock are shown in Figure 5. 
Values of k used in the JST and Turkel’s matrix 


dissipation model are given in §4.3. We used higher 
values of stencil Reynolds number for the Tam dis- 
sipation model to avoid excessive smearing of the 
shock profile. Both MacCormack-Baldwin (MB) 
and Turkel’s matrix dissipation model gave sharp 
shock profiles. Other models also performed rea- 
sonably well. The Tam dissipation model appears 
to be slightly better in eliminating preshock oscil- 
lations, but also seems to smear the shock profile 
more than other models. 

Our second test problem has a very sharp gra- 
dient in the initial condition. In Figure 6, we show 
the initial and computed profiles of this test prob- 
lem at time equal to 50. Computed solutions were 
obtained using the DRP and the standard 6th or- 
der central difference schemes. In this plotting scale, 
differences between the DRP and the 6th order cen- 
tral difference scheme are not noticeable. The Tam 
dissipation model with the stencil Reynolds number 
equal to 0.1 was used for these computations. There 
are minor oscillations near the shock. We concen- 
trate on that region in the next two figures. In Fig- 
ure 7, we compare solutions with the DRP scheme 
with and without artificial dissipation. Model con- 
stants k( 2 ) and in the matrix dissipation were 
kept as in the case shown in Figure 2. Both the Tam 
dissipation and Turkel’s matrix dissipation models 
were effective in eliminating spurious oscillation in 
the computed solutions. The Tam dissipation model 
with the stencil Reynolds number (R$) equal to 
0.05 appears to be slightly better that the recom- 
mended value i.e., 0.1 in removing most of the os- 
cillations. The shock location is slightly different 
with these two values of the stencil Reynolds num- 
bers. We then compare the DRP and the standard 
central difference schemes with the Tam dissipation 
model in Figure 8. Differences in the computed so- 
lutions with the DRP and the standard central dif- 
ference schemes appear to be small. This may be 
partially caused by the presence of large amount 
of artificial dissipation in the computed solutions. 
Next we examine the dissipation models in the Mac- 
Cormack scheme for computation of the test prob- 
lem 2 at time equal to 50. In Figure 9, we show 
the density profile computed using the MacCormack 
scheme without dissipation and compare it with the 
solution using the DRP scheme with the Tam dissi- 
pation model ( R g = 0.1). The MacCormack scheme 
without dissipation generates oscillations near the 
shock. We examined the effect of various dissipation 
models on the MacCormack scheme in Figure 10. 
The shock profile is slightly sharper with the ma- 
trix dissipation model than other models. As in the 
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case with the test problem 1, preshock oscillations 
were slightly less with the Tam dissipation model, 
but the shock profile was slightly more smeared. 
We examine the profile near the contact disconti- 
nuity in Figure 11. For computations in this re- 
gion, MacCormack- Baldwin (MB) and Turkel’s ma- 
trix dissipation model performed well. 

Our last test problem is a shock tube problem 
and was used by Lax^. In this paper we show 
computed solutions of this problem at time equal 
to 0.16. For this model problem, the values of k 
used in the matrix dissipation model are 16 times 
the recommended values of the JST model (§4.3). 
This choice was made only to ensure sufficient arti- 
ficial dissipation for the test cases without any op- 
timization of the values of k. As shown in Figure 
12, major differences in computed solutions are lim- 
ited to regions around the discontinuities. We fo- 
cus on this region for computations with the DRP 
and MacCormack schemes in next two figures. In 
Figure 13 we compare computed density profiles us- 
ing the DRP scheme with the exact solution. The 
computed profile with the stencil Reynolds num- 
ber equal to 0.1 appear to be in reasonable agree- 
ment with the exact solution. Solution using the 
fourth order DRP scheme (seven point stencil) is 
very close to that using the standard sixth order 
central scheme (seven point stencil). For the sten- 
cil Reynolds number equal to 0.05, the solution is 
smeared and the shock location is slightly shifted. 
The smearing of the shock profile is due to excessive 
amount of artificial disspation. The change in the 
shock location may in part be caused by the Tam 
dissipation model which is not cast in conservation 
form. Tam and Shen^, however, pointed out that 
some dispersion error may come from their damping 
model. We next examine computations using the 
fourth order MacCormack Scheme. Solutions using 
this scheme are shown in Figure 14. As in the case 
of previous test problems, computed solutions with- 
out any artificial dissipation exhibit oscillations. We 
also show the results with the MacCormack- Baldwin 
and Turkel’s matrix disspation model in this figure. 
Artificial dissipation models eliminated most of the 
spurious numerical oscillations. Our choices of the 
values of k may have caused excessive damping with 
the matrix dissipation model. Finally we show com- 
puted pressure profiles of Lax’s problem in Figures 
15 and 16. Here again we find that artificial dis- 
sipations eliminated most of the oscillations in the 
DRP (see figures 2, 7 and 8 for solutions using the 
DRP scheme without artificial dissipation) and the 


fourth order MacCormack scheme without artificial 
dissipation. Basic observations are essentially the 
same as from computed density profiles. 

In summary, we would like to say that various 
artificial models considered in this study in general 
performed well. However, one probably should opti- 
mize some parameters to improve the quality of the 
solutions. The fourth order MacCormack scheme 
without any artificial dissipation generates spuri- 
ous oscillations near sharp gradients. We recom- 
mend adding artificial dissipation for computations 
of nonlinear waves or flows with shocks with this 
scheme. Artificial dissipation plays a more impor- 
tant role in getting good computed solutions using 
the DRP scheme. Without any artificial damping, 
the DRP scheme appears to generate excessive nu- 
merical oscillations. 
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Figure Is Evolution of the density pulse in problem 1 



Figure 2: Effect of dissipation on the DRP scheme for problem 1 
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Figure 3: Comparison of central schemes using the Tam dissipation model for problem 1 
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Figure 5: Effect of artificial dissipation on the MacCormack scheme for problem 1 









X 

Figure 7: Effect of artificial dissipation on the DRP scheme for problem 2 
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Figure 8: Comparison of central schemes using the Tam dissipation model for problem 2 
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Figure 15: Computed pressure profile of Lax’s problem 
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Figure 16: Computed pressure profile near the shock 
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